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Production of muons and neutrinos in cosmic ray interactions with the atmosphere has been investigated with 
a cascade simulation program based on Lund Monte Carlo programs. The resulting 'conventional' muon and 
neutrino fluxes (from 7r, K decays) agree well with earlier calculations, whereas the improved charm particle 
treatment used in this study gives significantly lower 'prompt' fluxes compared to earlier estimates. This implies 
better prospects for detecting very high energy neutrinos from cosmic sources. 



1. Introduction 

The flux of muons and neutrinos at the earth 
has an important contribution from decays of par- 
ticles produced through the interaction of cosmic 
ray particles in the atmosphere. This has an in- 
terest in its own right, since it reflects primary in- 
teractions at energies that can by far exceed the 
highest available accelerator energies. It is also 
a background in studies of neutrinos from cos- 
mic sources as attempted in large neutrino tele- 
scopes, such as Amanda, Baikal, Dumand and 
Nestor. Here we give a short report of a compre- 
hensive study of muon and neutrino production 
in cosmic ray proton interactions with nuclei in 
the atmosphere using detailed Monte Carlo sim- 
ulations [1,2]. 

The atmospheric muons and neutrino fluxes are 
dominated by 'conventional' sources, i.e. decays 
of relatively long-lived particles such as 7r and 
K mesons. This is well understood from ear- 
lier studies [3,4] and confirmed by our investiga- 
tions [1,2]. With increasing energy, the proba- 
bility increases that such particles interact before 
decaying. This implies that even a small frac- 
tion of short-lived particles can give the dominant 
contribution to the high energy tail of the muon 
and neutrino fluxes. These 'prompt' muons and 
neutrinos arise through semi-leptonic decays of 
hadrons containing heavy quarks, most notably 
charm. 
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Previous estimates of the flux of prompt muons 
and neutrinos from charm [5-7] vary by a few 
orders of magnitude (cf. Fig. 3). This is due to 
the different models used to calculate the charm 
hadron cross section and energy spectra. The 
models are usually constrained or fitted to ex- 
perimental data at accelerator energies and then 
extrapolated to the higher energies needed in this 
context. Obviously, the extrapolation can only be 
trustworthy if based on a sound physical model. 
The main new contribution of our study is in this 
context, where we apply state-of-the-art models 
to simulate charm particle production in high en- 
ergy hadron-hadron interactions. 

2. The cascade algorithm 

The prompt flux is essentially independent of 
the zenith angle due to the negligible probability 
that the very short-lived particles interact before 
decaying. This flux is therefore relatively more 
important in the vertical direction where the con- 
ventional flux is lower. The fluxes are therefore 
calculated for the vertical direction which has the 
advantage that a simple exponential model for the 
density of the atmosphere can be used, i.e. the 
atmospheric depth is 

X[g/cm 2 ] = 1300 e"^ 1 "]/ 6 - 4 (1) 

The parameters are chosen to give a good descrip- 
tion of the atmosphere at altitudes of 10-40 km 
where the interactions typically take place. 



The simulation algorithm is constructed as fol- 
lows [1,2]. The energy of a cosmic ray proton 
is chosen from the form Q(E) ~ E~ 1_1 , with 
7= 1.7 for E below 5-10 6 GeV and 7 = 2 above. A 
proton-nucleon interaction is then generated with 
Pythia [8] resulting in a complete final state of 
particles. Cross section and mean free path are 
scaled according to the composition of the at- 
mosphere, and an interaction height is chosen. 
The produced particles are traced through the 
atmosphere until they either decay or interact 
as given by simulated decay lengths and interac- 
tion lengths. Particle decays are fully simulated 
with daughter particle momenta. In case of in- 
teractions, the interacting particle is regenerated 
in the same direction but with degraded energy 
corresponding to the appropriate leading parti- 
cle spectrum. The procedure is repeated until all 
particles have decayed, hit the ground or their en- 
ergy fallen below a minimum of 10 2 GeV. Finally, 
muons and neutrinos are counted with appropri- 
ate weights to obtain their energy spectra. 

The simplification that secondary interactions 
are not simulated in full detail is justified since 
only leading particles will contribute significantly. 
Other particles with lower energy give a negligi- 
ble flux compared to the flux of primary parti- 
cles at the corresponding energy. Furthermore, 
proton-nucleus and meson-nucleus collisions are 
very similar except for the leading particle. 

To obtain the flux of conventional muons and 
neutrinos, an inclusive event sample is simulated 
using Pythia in a mode generating minimum 
bias proton-proton interactions (including diffrac- 
tive scattering). The particle production results 
from Lund model hadronisation of colour string 
fields between partons scattered in semi-soft QCD 
interactions. The prompt flux is obtained from 
charm particles arising from the hadronisation of 
charm quarks produced in the processes gg — > cc 
and qq — > cc as calculated with leading order per- 
turbative QCD (pQCD) matrix elements. 

3. Numerical results 

From the extensive studies in [1,2] we show 
in Fig. 1 the main results in terms of the calcu- 
lated fluxes of muons and neutrinos. The prompt 
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Figure 1. Vertical fluxes (_E 3 -weighted energy 
spectra) of muons and neutrinos from conven- 
tional (7r, K decays) and prompt (charm decays) 
sources and their sum as obtained from simulated 
cosmic ray interactions in the atmosphere. 



contribution from charmed particles dominates at 
high energies. This contribution is subdivided on 
the different charmed particles in Fig. 2, demon- 
strating the dominance of the D^' mesons. 

The results for the inclusive prompt and con- 
ventional fluxes can be parametrised as (cf. the 
primary flux) 

( NqE-i- 1 + AE) E < E 
4>(E)={ * ' (2) 

[ Nl>E-i'- 1 /{l + A' E) E>E 



Table 1 



Values of parameters in eq. 2 obtained from fits to the Monte Carlo results in Fig. 1. 









N 




7 




A 




Eq 


i 




A 1 




K 




Conv. /i~ 


+ M + 


1.3 


■ io- 


1 


1.77 


4.4 


■ 10" 


3 


2.7- 10 5 


2.40 


9.2 


■ io- 


6 


4.6 • 10- 


1 


Prompt /i~ 


+ M + 


2.4 


• io- 


6 


1.63 


6.5 


• 10" 


7 


3.2 • 10 7 


2.64 


6.2 


• io- 


8 


1.2 • 10 1 




Conv. - 


f f/i 


8.5 


• 10" 


3 


1.55 


9.1 


• io- 


3 


7.9 • 10 3 


2.69 


1.0 


• 10" 


6 


3.3 • 10- 


4 


Prompt - 




2.5 


• io- 


6 


1.63 


6.5 


• 10" 


7 


4.4- 10 7 


2.78 


1.3 


• 10" 


8 


8.0 • 10 1 




Conv. v e - 


- v e 


3.8 


• 10" 


4 


1.70 


3.9 


• 10" 


3 


1.3 • 10 5 


2.32 


4.6 


• 10" 


6 


1.8 • 10- 


3 


Prompt v e - 


- v e 


2.5 


• 10" 


6 


1.63 


6.5 


• io- 


7 


4.1 • 10 7 


2.87 


6.0 


• 10" 


7 


6.4 • 10 3 






4 5 
log 10 (E/GeV) 

Figure 2. Flux of muon neutrinos and anti- 
neutrinos from decays of the different charmed 
particles. 



with an accuracy of typically better than 10% 
using the fitted parameter values in Table 1. 

4. Discussion 

Whereas earlier calculations [3,4] agree well 
with our conventional muon and neutrino fluxes 
[1,2], previous estimates [6,7] of the prompt fluxes 
are up to orders of magnitude larger as illustrated 
in Fig. 3. (Prompt fluxes are direction indepen- 
dent up to ~ 10 7 GeV [1] and therefore directly 
comparable between horizontal and vertical di- 
rections.) These differences are due to different 
models for charm production, both regarding the 
energy dependence of the cross section and the 



distribution in longitudinal momentum fraction 
xp = p\\/pmax (Feynman-a;) of the charmed par- 
ticles. Model I [6] uses a parametrised energy de- 
pendence of the cross section, normalised to some 
charm production data, and an x j?-spectrum of 
the form 



dN/dx F ~ (i - x F y 



(3) 



with ajj = 5 and «a c = 0.4 for _D-mesons and A c . 
Feynman scaling is thus assumed, i.e. xp is in- 
dependent of the collision energy. Model II [7] 
simply takes charm as 10% of the total inelastic 
cross section and xp distributions of eq. (3) with 
ajj = 3 and «a c = 1- Model III [7] corresponds 
to charm quark production calculated with lead- 
ing order pQCD matrix elements using relatively 
hard parton distributions. 
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Figure 3. Prompt and conventional muon fluxes 
(solid lines) compared with earlier model calcula- 
tions I— III of horizontal prompt fluxes (see text). 




Figure 4. Distribution in xp of _D-mesons in 
a model obeying Feynman scaling (ajj = 5 in 
eq. (3), solid line), and from Pythia at pro- 
ton beam energy 10 3 GeV (dashed line), 10 6 GeV 
(dotted line) and 10 9 GeV (dash-dotted line). 



Whereas our model reproduces available data 
on charm production cross sections [2], not all of 
models I— III do. When extrapolating orders of 
magnitude in energy, the cross section differences 
are enhanced to factors of 10. A more important 
reason for our lower flux is, however, the strong 
breaking of Feynman scaling as demonstrated in 
Fig. 4. This arises from the dominance of charm 
production close to threshold in pQCD such that, 
with increasing energy, smaller momentum frac- 
tions x of incoming partons will contribute (giving 
a sensitivity to parton density parametrisations). 
Together with the hadronisation mechanism, this 
dynamics imply smaller x p at higher collision en- 
ergies. 

Finally, the calculated atmospheric muon- 
neutrino flux is considered as a background in 
comparison to astrophysical sources, as illus- 
trated in Fig. 5. The neutrino production in cos- 
mic ray interactions with the galactic medium 
and, in particular, the diffuse neutrino flux from 
active galactic nuclei is estimated to be in excess 
of the atmospheric neutrino background for high 
energies and should therefore be possible to de- 
tect with large scale neutrino telescopes. 
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Figure 5. Vertical fluxes of conventional and 
prompt atmospheric muon-neutrinos (solid lines) 
compared to astrophysical sources: the flux from 
cosmic ray interactions with the intra-galactic 
medium (dashed curves derived from [9]; upper 
curve apply for the direction towards the center 
of the galaxy, and lower curve orthogonal to the 
galactic plane) and the estimated diffuse fluxes 
from active galactic nuclei (dotted line [10], dash- 
dotted line [11]). 
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